!! Existing code:
!!
![readinfo]
!  [readgrids]
!    read (infile, *) nz
!    allocate (zj(nz))
!    read (infile, *) zj(1:nz)
!  [read_ecpars]
!    allocate (piz(nz,nz))
!    do jz=1,nz
!       read (infile, *) (piz(jz,jzpr), jzpr=1,nz)  !transition
!    end do
!! *****************************************************
!! Calling from code in readinfo
!
!use tauchen_mod
! ! .. discretize persistent inc process z ..
! call tauchen(zero, sd_eps, rho_z, multiple_sd, nz, zj, piz)
!
!! Need specify new parameters to be read in [read_ecpars]
!sd_eps
!rho_z
!multiple_sd
!
!! Move readinf nz to [read_ecpars]
!

! ********************************************
module tauchen_mod

implicit none

integer, parameter ::  ik=selected_int_kind(9)
integer, parameter ::  rk=kind(1.d0) !!!! prec

private
public:: tauchen, simulate_markov !  linspace

contains
!   [tauchen]
!       call linspace
!       normal
!   <normal>
!       call normaldensity
!       contains
!           [normaldensity]


! G. Tauchen (1986) 'Finite State Markov-Chain Approximations to
! Univariate and Vector Autoregressions' Economics Leters 20: 177-181

! Contains numerical integration of normal density to yield distribution.
subroutine tauchen(meaninnov, stdinnov, persistence, multiple, znum, zvec, pi)

integer(ik):: znum, j, gridsize, k
real(rk)::meaninnov, stdinnov, persistence, multiple, zvec(znum), pi(znum, znum), &
		  f1, f0, stdz, zlow, zhigh, meanz, w, z, lowerbound

intent(in):: meaninnov, stdinnov, persistence, multiple, znum
intent(out):: zvec, pi

! Set endpoints of stochastic variable grid at multiple m of the
! standard deviation.

stdz = stdinnov**2.0_rk
stdz = stdz/(1.0_rk - persistence**2.0_rk)
stdz = dsqrt(stdz)
meanz = meaninnov/(1.0_rk - persistence)
zlow = meanz - stdz*multiple
zhigh = meanz + stdz*multiple

lowerbound = meaninnov - stdinnov*dmax1(10.0_rk, 2.0_rk*multiple)
gridsize = 10000

call linspace(zlow, zhigh, znum, zvec)

pi = 0.0_rk

w = (zhigh - zlow)/dble(znum-1_ik)

! this is equations (3a) and (3b) from tauchen (1986)

do j = 1, znum
	
	z = zvec(1) - persistence*zvec(j)
	f1 = normal(z + w/2.0_rk, meaninnov, stdinnov, gridsize, lowerbound)
	pi(j,1) = f1
	
	do k = 2, znum - 1
		z = zvec(k) - persistence*zvec(j)
		f1 = normal(z + w/2.0_rk, meaninnov, stdinnov, gridsize, lowerbound)
		f0 = normal(z - w/2.0_rk, meaninnov, stdinnov, gridsize, lowerbound)
		pi(j,k) = f1 - f0
	end do

	z = zvec(znum) - persistence*zvec(j)
	f0 = normal(z - w/2.0_rk, meaninnov, stdinnov, gridsize, lowerbound)
	pi(j,znum) = 1.0_rk - f0

end do

end subroutine tauchen

function normal(upperbound, mean, sd, gridsize, lowerbound)

! Program uses numerical integration with gridsize grid points to
! compute the pr(lowerbound < x < upperbound) for a normal distri-
! bution with mean given and standard deviation of sd.  Numerical
! integration is based on a simple upper and lower darboux sum
! method where the integrand is the normal density function.

! Aubhik, bastion of stupidity, 18.04.2002

optional:: mean, sd, gridsize, lowerbound
integer(ik):: gridsize
real(rk):: lowerbound, upperbound, mean, sd, increment, normal
real(rk), allocatable:: x0dl(:), x1du(:), f0dl(:), f1du(:), f(:)

! Arguments for mean, sd, grid size for numerical integration and the
! lower bound of the interval are optional.
if (present(mean)) then
	continue
else							! set mean to 0 if not specified	
	mean = 0.0_rk
end if

if(present(sd)) then
	continue
else							! set standard deviation to 1 is not specified
	sd = 1.0_rk
end if

if(present(lowerbound)) then
	continue
else							! give probability from roughly - infinity	
	lowerbound = -10.0*sd + mean
end if


if(present(gridsize)) then
	continue					! use a 1000 points for numerical integration if not speficied
else
	gridsize = 100000_ik			
end if

! Grids for lower and upper darboux sums, remember Megan and math analysis at DRL.
increment = (upperbound - lowerbound)/dble(gridsize)

allocate(x0dl(gridsize), x1du(gridsize), f0dl(gridsize), f1du(gridsize), f(gridsize))

call linspace(lowerbound, upperbound - increment, gridsize, x0dl)
call linspace(lowerbound + increment, upperbound, gridsize, x1du)

! The normaldensity subroutine gives the value of the normal density at each grid point.
call normaldensity(gridsize, x0dl, mean, sd, f0dl)
call normaldensity(gridsize, x1du, mean, sd, f1du)

! Average the function values.
f = (f0dl + f1du)/2.0_rk

! Weight each one by the width of the rectangle of integration.
normal = sum(f*increment)

contains

subroutine normaldensity(numberofpoints, vectorofpoints, mean, sd, densities)

integer(ik):: numberofpoints
real(rk):: vectorofpoints(numberofpoints), densities(numberofpoints), mean, sd, &
		   variance, pi, coefficient, transcend(numberofpoints)

intent(in):: numberofpoints, vectorofpoints, mean, sd
intent(out):: densities

! From the matlab version of this program, the reference is apparently
! page 244, chapter 4 of bopwerman and o'connell (1997) applied statistics

variance = sd**2.0_rk
pi = 2.0_rk*dasin(1.0_rk)
coefficient = variance*2.0_rk*pi
coefficient = 1.0_rk/dsqrt(coefficient)

! Note the vector-valued operations as a scalar, mean, is subtracted from each
! element of a vector, then this resultant vector is divided by two scalars and
! finally used as a vector of exponents for the exponential function.

transcend = (vectorofpoints - mean)**2.0_rk
transcend = transcend/variance
transcend = transcend/2.0_rk
transcend = -1.0_rk*transcend
densities = coefficient*dexp(transcend)

end subroutine normaldensity

end function normal


!!!!!!!!!!!!!!!!!
!	LinSpace	!
!!!!!!!!!!!!!!!!!

subroutine linspace(lb, ub, gridnum, x)

integer(ik):: gridnum, j1
real(rk):: lb, ub, x(gridnum), y(gridnum-2_ik)

intent(in):: lb, ub, gridnum
intent(out):: x

! This subroutine is written to reproduce the MATLAB 4.2 func-
! tion linspace.  It generates a vector which contains a grid
! of n equally spaced elements between loweround and upper-
! bound.  Note this version produces a row vector.

do j1 = 1_ik, gridnum - 2_ik
	y(j1) = dble(j1)
end do

! Note that I am multiplying an n-2x1 vector, Y, by a scalar
! (which implies multiplying each element of the vector by
! this scalar) and then adding another scalar (which implies
! adding this scalar to each element of the vector).

y = ((ub - lb)/(gridnum-1_ik))*y + lb

x(1_ik) = lb
x(2_ik:gridnum-1_ik) = dble(y)
x(gridnum) = ub

end subroutine linspace


! *******************************************************
!
! SIMULATE MARKOV CHAIN
!
! *******************************************************

subroutine simulate_markov(gamma,logzj,nz,ntime)

 implicit none
 integer :: j, i

 integer :: nz
! real(rk), dimension(:), allocatable:: logzj(:)
! real(rk), dimension(:), allocatable:: gamma(:,:)

 real(rk), dimension(:):: logzj(:)
 real(rk), dimension(:):: gamma(:,:)

 integer, parameter:: nshocks=1
 integer :: izpr, iz, time, ntime

 real(rk), dimension(:), allocatable:: ran_out(:)
 integer, dimension(:), allocatable:: seed(:)
 real(rk), dimension(:), allocatable:: logz_t(:)

 real(rk) :: parsum
 real(rk) :: x_z
 real(rk) :: logz_markov_aver, logz_markov_autocorr, logz_markov_sd


! allocate (logzj(nz))
! allocate (gamma(nz,nz))
 allocate (ran_out(nshocks))
 allocate (seed(nshocks))
 allocate (logz_t(1:ntime))



! Initial States
    iz=1
 seed = (/ 7622321/)
    CALL random_seed(put=seed(1:nshocks))

 do time=1, ntime
    CALL random_number(ran_out(1:nshocks))
    x_z = ran_out(1)
    ! .. Update income states ..
    if(x_z.eq.1._rk) then
        izpr=nz
    else if(x_z.eq.0.0_rk) then
        izpr=1
    else
      parsum=0.0_rk
      do izpr=1, nz
       if((x_z.ge.parsum).and.(x_z.lt.parsum+gamma(iz,izpr)))  exit
       parsum=parsum+gamma(iz,izpr)
      end do
    end if
    logz_t(time)=logzj(iz)
    !  update
     iz=izpr
 end do !time


!    open(297, file='markov_data.dat')
!     do time=1, ntime/10
!        write(297,'(i6,f14.6)') time, logz_t(time)
!     end do
!    close(297)

    open(297, file='markovsimul.dat')
    write(297,*)
    write(297,'(a30)') 'STATISTICS MARKOV CHAIN '
    logz_markov_aver=sum((logz_t(1:ntime)))/real(ntime)
    write(297,'(a18,f9.6)') 'mean logz=', logz_markov_aver
    logz_markov_sd=sqrt(sum(((logz_t(1:ntime))-logz_markov_aver)*&
           & ((logz_t(1:ntime))-logz_markov_aver))/real(ntime))
    write(297,'(a18,f9.6)') 'sd logz=', logz_markov_sd
    logz_markov_autocorr= &
     & (sum(((logz_t(2:ntime))-logz_markov_aver)* &
        &   ((logz_t(1:ntime-1))-logz_markov_aver))/real(ntime-1) &
        & )/ &
        & (logz_markov_sd*logz_markov_sd)
    write(297,'(a18,f9.6)') 'autocorr logz=', logz_markov_autocorr
    close(297)

 deallocate(logz_t)

end subroutine simulate_markov








end module tauchen_mod
